Density profiles of loose and collapsed cohesive granular structures generated by 

ballistic deposition 
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Loose granular structures stabilized against gravity by an effective cohesive force are investigated 
on a microscopic basis using contact dynamics. We study the influence of the granular Bond number 
on the density profiles and the generation process of packings, generated by ballistic deposition under 
gravity. The internal compaction occurs discontinuously in small avalanches and we study their size 
distribution. We also develop a model explaining the final density profiles based on insight about 
the collapse of a packing under changes of the Bond number. 

PACS numbers: 47.57.-s, 45.70.Mg, 83.80.Hj 



I. INTRODUCTION 



Loose granular packings, met ast able granular struc- 
tures and fragile granular networks play an important 
role in a wide range of scientific disciplines, such as 
collapsing soils [H-Q^ ^ ne powders [E| or complex flu- 
ids p, 0- In collapsing soils without any doubt there 
is a metastable or fragile granular network involved [l]— 
H, H, 111 • A similar failure behavior can be found in col- 
loidal gels [lo|| and snow [ll|, [TgJ . But also powders have 
in most cases an effective cohesive force, e.g. due to a cap- 
illary bridge between the particles or van der Waals forces 
(important when going to very small grains, e.g. nano- 
particles) leading to the formation of loose and fragile 
granular packings H, 0, [l3| . In many complex fluids a 
fragile/metastable network of colloids/grains is believed 
to be the essential ingredient for the occurrence of shear 
thickening \^ or yield stress behavior Q. 

The general feature of such fragile networks is that they 
can collapse/ compact under the effect of an applied load 
0, EL • This load can be an external load or exerted 
internally by a force acting on all particles within the 
structure. This "internal collapse" is important in dif- 
ferent applications like cake formation of filter deposits 
[l4]-[l6| , where the compaction force in most situations is 
the drag force exerted on the grains by the flow which 
is typically porosity dependent [17]. The structure's own 
weight leads to compaction of snow after deposition [l8| 
and during aging [l9|, [20j, or to sediment compaction 
[2ll-[23|. In all cases, typically a depth dependent poros- 
ity is observed and quantified by continuum descriptions 
US Ql| l2ll - [23| . In most cases the details of the poros- 
ity profile are influenced by a combination of different me- 
chanical and chemical processes [HI, [2l| . It is well known 
that the porosity of a structure is of major importance 
for its mechanical properties [ll|, [2l|, Ell ? m filtration 
processe s I26j and its chemical properties like catalytic 
activity [27| • The aim of this paper is to study the micro- 
scopic processes, i.e. on the grain scale, for these internal 
compaction processes. For this, we will investigate the 
compaction due to gravity in a simplified model system of 
grains held together by cohesive bonds. We analyze how 



the density profiles depend on the granular Bond num- 
ber, i.e. the ratio of cohesive force by gravity, and what 
the influence of the dynamics of deposition/collapse. As 
discussed above loose structures are generated in nature, 
industrial application, experiments or simulation by dif- 
ferent processes. Here, we focus on ballistic deposition. 
However, we expect the findings of this paper to be of 
relevance to all systems involving compaction due to the 
particles' own weight. 

After a description of the simulation model and a brief 
discussion of possible experimental realizations in sec- 
tion [III we first study the resulting density profiles when 
gravity acts during deposition, in particular the influ- 
ence of the granular Bond number (sec. 1111)) . To under- 
stand the shape of the density profiles we study in the 
following (sec. IIV|) the role of the dynamics of the col- 
lapse occurring in small avalanches. We study the aver- 
age "avalanche profile" defined here as the average dis- 
tance a particle moves downwards after being deposited 
depending on its height. We observe characteristic pro- 
files which can be used to relate the final density pro- 
file to the deposition density, given by the number of 
deposited particles per unit volume (sec. [V)). To under- 
stand this phenomenological profile we study a simpler 
system where first all particles are deposited followed by 
the collapse of the whole structure leading to an even 
simpler profile (sec. IVI|) . Our calculations yield that this 
linear profile is obtained in all processes where a homo- 
geneous initial configuration is collapsed/compacted to a 
homogeneous final state. In Sec JVIII we show that the 
phenomenological obtained avalanche profile obtained in 
sec. IIIII can be derived from the linear avalanche profiles 
of the homogeneous collapse. 



II. DESCRIPTION OF SIMULATION MODEL 

The dynamical behavior of the system during gener- 
ation is modeled with a particle based method. Here 
we use a two dimensional variant of contact dynamics, 
originally developed to model compact and dry systems 
with lasting contacts [28l43lj. The absence of cohesion 
between particles can only be justified in dry systems 
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on scales where the cohesive force is weak compared to 
the gravitational force on the particle, i.e. for dry sand 
and coarser materials, which can lead to densities close 
to that of random dense packings. However, an attrac- 
tive force plays an important role in the stabilization of 
large voids [32j, leading to highly porous systems as e.g. 
in fine cohesive powders, in particular when going to very 
small grain diameters. Also for contact dynamics a few 
simple models for cohesive particles are established [32I — 
|35| ]. Here we consider the bonding between two particles 
in terms of a cohesion model with a constant attractive 
force F c acting within a finite range d c , so that for the 
opening of a contact a finite energy barrier F c d c must 
be overcome. In addition, we implement Coulomb and 
rolling friction between two particles in contact, so that 
large pores can be stable [32|, l36l - l39| . 

To generate the loose structure we use ballistic depo- 
sition where each deposited particle, chosen at random 
horizontal position, is attached to the structure at maxi- 
mal possible height with zero velocity. At the same time 
we allow for all particles to move which can lead to a par- 
tial collapse of the structures due to gravity [M@,fl^l40|. 
The structure is deposited on a fiat surface, i.e. a wall at 
the bottom. We use periodic boundaries in horizontal 
direction to avoid effects of side walls, like Janssen effect. 
During this process the time interval between successive 
depositions crucially determines the structure and den- 
sity profiles of the final configurations. Here we will focus 
on the two extreme cases of very large time intervals, i.e. 
the system can fully relax after each deposition of a single 
grain, and vanishing time interval, i.e. the collapse of the 
systems happens after the deposition process is complete. 
In the first case the interval is chosen large enough to let 
the system compactify and relax due to the additional 
weight of the deposited grain. This is verified on the one 
hand by checking that the final density is independent on 
the time interval and on the other hand by monitoring 
the dynamics of the process. Having no time between 
depositions in practice means that first we perform pure 
ballistic deposition [4l|, |42|, and then switching on the 
full particle dynamics leading to a collapse of the system 
due to gravity. Experimentall y, t he two cases can be re- 
alized in a Hele-Shaw cell |43i-l45j which can be tilted to 
effectively change gravity In the slow deposition process, 
simply the cell is slowly filled in an upright position so 
that full gravity acts on the grains. In the other case 
the Hele-Shaw cell will be almost horizontal, so that the 
grains can be filled in with nearly vanishing gravity, and 
then the cell is tilted so that gravity can fully act on the 
grains, leading to an abrupt collapse of the structure. 



III. DENSITY PROFILES WHEN GRAVITY 
ACTS DURING DEPOSITION 

In this section we analyze the density profiles for the 
case of large enough time intervals between successive 
depositions to allow the systems to relax under the ef- 



fect of gravity as described in the previous section. It is 
expected that the density and the characteristics of the 
density profiles are mainly determined by the ratio of the 
cohesive force F c to gravity F g , typic ally defined as the 
granular Bond number Bo g = F c /F g [46|, E3- Obviously 
the case of Bo g = corresponds to the cohesionless case 
whereas for Bo g — >> 00 gravity is negligible. A similar 
dimensionless quantity had been identified as most im- 
portant parameter in previous studies on compaction of 
cohesive powders [32|, S, |48j . 

In the following, we use monodisperse systems with a 
friction coefficient \i = 0.3 and a rolling friction coeffi- 
cient of jjL r = 0.1 (in units of particle radii). The effect of 
varying these parameters is also studied exemplary and 
will be discussed later. Typically the values of the den- 
sity can depend on these parameters as shown in Ref. 
[37j whereas the qualitative behavior does not change. 
Figure Q] shows the final structures obtained for differ- 
ent values of granular Bond number ranging from to 
10 6 . Also the limit of infinite Bond number is shown, 
leading to pure ballistic deposition [42j well studied al- 
ready in the past. For small Bond numbers, here repre- 
sented by Bo g = 0, the system typically reaches a ran- 
dom close packing which also has been studied intensively 
in the past. Note that our case of monodisperse parti- 
cles typically leads in dense packings to crystallization 
effects which could be avoided by using a small poly- 
dispersity. As our focus in this paper is on the looser 
structures where this effect is not very important we pre- 
fer the monodisperse system to keep the model as simple 
as possible. In the intermediate range of Bond numbers 
the density varies between the two limiting values. 

Plotting the density profile depending on the vertical 
position y (Fig. [2]) provides a more quantitative analysis. 
It can be seen that in the two limiting cases (Bo g = and 
Bo g — >> 00) the density is constant. For the infinite Bond 
number this can be explained easily as no collapse at all 
occurs and the density profile is that of a ballistic depo- 
sition and thus constant [4l|, |42|. For the non-cohesive 
case a close packing is expected, also leading to a con- 
stant density. This will be discussed again in more detail 
later in this paper. In the intermediate range the density 
decreases with increasing height. This is a result of the 
generation process where the fragile structure is partially 
collapsed due to the weight of the added particles which 
happens discontinuously in relatively small avalanches as 
will be discussed in more detail in the next section (sec. 

IYJ). 

Knowing that the density depends on vertical position 
a general dependence of the total density on the Bond 
number cannot easily be defined. Instead, for a given 
system size as in Fig. [2] the density at a fixed position 
can be measured. In Fig. [3] the averaged density in the 
lower half excluding the region very close to the bottom 
is shown versus the granular Bond number. The den- 
sity varies between the two limiting cases Bo g = and 
Bo g — >• 00. Note that the Bond number is plotted in a 
logarithmic scale, i.e. to see substantial changes of vol- 
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FIG. 1: (Color online) Final structures achieved by the deposition/collapse process for different granular Bond numbers Bo g . 
In addition to the particles compressive forces are illustrated by red (dark gray) lines connecting the center of masses between 
the particles. In the case of Bo g — >• oo no forces are present as it is realized in the simulations by switching off gravity. 
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FIG. 2: (Color online) Density profiles for different granular 
Bond numbers Bo g (cf. fig. [lj - Here, the volume fraction v 
is plotted. In this case the volume fraction is measured in 
thin slices of given width (here: 3.97 particle radii) at varying 
height y. For Bo g — no cohesion is active and the random 
close packing is reached. In the limit Bo g — » oo the system 
does not collapse at all, and the simple ballistic deposition 
case [43] is obtained. 



FIG. 3: (Color online) Average volume fraction Slower half 
depending on granular Bond number Bo g . The density is 
averaged in the lower half of the system excluding the region 
very close to the bottom to avoid border effects (here we 
excluded the region below the height of 50 particle radii so 
that clearly boundary effects are removed for all curves, cf. 
fig. [2J. The volume fractions vary between the two limits 
given by random close packing (Bo g = 0) and pure ballistic 
deposition (Bo g — > 00). 



ume fraction the cohesive force or the gravitational force 
have to be changed by orders of magnitude. Particles 
with similar gravity and cohesive force will show the same 
typical behavior. As typically both forces depend on the 
size of the particles it appears to be natural to character- 
ize the behavior of granular matter and powders by the 
grain size. For non-cohesive material recent experimen- 
tal, numerical and theoretical studies [49l-[52j investigate 
the influence of the friction coefficient on, e.g. the volume 
fraction. A similar behavior as found here for the co- 
hesive material when varying the granular bond number, 



has been found [HI, [50[ : varying the friction coefficient 
on a logarithmic scale leads to a variation between the 
values 0.84 for the packing fraction of a random close 
packing and the value 0.77 for infinitely large friction co- 
efficient (in two dimensions, in three dimensions between 
0.64 and 0.55). In the cohesive case as discussed here 
this range of accessible volume fractions is much higher 
and limited by the preparation protocol, i.e. in this paper 
by the ballistic deposition. This limit of course can be 
changed when changing the preparation protocol, e.g. by 



introducing a capture radius (cf. sec. PVT]) . 
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FIG. 4: (Color online) Illustration of the effect of system size 
for intermediate density range (Bond number Bo g — 10 3 ). 
Plotting the depth H — y measured from the surface H of the 
final packings smaller systems show the same profile as large 
systems. 
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FIG. 5: Particle trajectories of particles of the small system 
(H = 147, cf. Fig. [4|) for the whole deposition/collapse. For 
better visibility only each 5th particle's trajectory is shown, 
i.e. the trajectories of 640 particles (instead of all 3200 par- 
ticles). Viewing the total system (and on the right more de- 
tailed when zooming in) illustrates that parts of the system 
move collectively downwards accompanied by a sidewards mo- 
tion or rotation. When zooming in the individual trajectories 
can be identified which are composed of the sum of paths 
during all the small avalanches experienced by the particle. 



For all results presented above the total system height 
H was fixed, i.e. the deposition process stops when no 
more particles can be deposited below a specified value 
H. When comparing density profiles for different system 
heights H plots depending on the vertical position y will 
show different densities. A scaling can be achieved when 
plotting the density versus the depth H — y as illustrated 
in Fig. HJ This means the upper part of the large system 
is depositing and collapsing in the same way as the small 
system while additionally leading to a further collapse of 
the structure deposited previously below, accompanied 
by a downwards motion of the whole upper part. Obvi- 
ously the slow deposition process guarantees that inertia 
is not important (cf. sec. PVT]) . 

The specific behavior of the density profiles shown in 
this section results from a deposition process combined 
with a collapse of the current structure due to gravity. 
The deposition is characterized by the number of de- 
posited particles per volume, which we call "deposition 
density" and which here is not constant (sec.[V)). The col- 
lapse happens successively in relatively small avalanches, 
analyzed in detail in the following section. In section IVl 
we will show that these avalanches can be used to relate 
the final density profile to the "deposition density" . 



IV. ANALYSIS OF THE AVALANCHES 
DURING DEPOSITION/COLLAPSING 

Typically the collapsing of the structures, as mentioned 
earlier, happens discontinuously in small avalanches. As 
these avalanches are important also for the final den- 
sity profiles (see sec. IIII|) their characteristics is stud- 
ied in detail in this section. To illustrate the nature of 
these avalanches in Fig. [5] the trajectories of the parti- 



cles are plotted for a relatively small system of height 
H = 147 consisting of about 3200 particles (for better 
visibility only each 5th trajectory is shown, i.e. the tra- 
jectories of 640 particles, instead of all 3200 particles). 
The avalanches are a collective motion of parts of the 
system. This mainly downwards motion is accompanied 
by a sidewards motion or rotation. When zooming in 
individual trajectories can be identified. These trajec- 
tories represent the motion of each particle during de- 
position/collapsing. Thus, they show the paths that a 
particle experiences in all avalanches at different times. 
Neighboring particles can have very similar trajectories, 
i.e. they belong to the same set of avalanches at different 
times. 

In Figure [6] we show the size of avalanches depending 
on initial and final vertical position. This size is mea- 
sured by Ay, the total downwards displacement of the 
particle after its deposition, i.e. initial position minus fi- 
nal position. This represents for each particle the sum of 
all avalanches occurring during the generation process, 
resulting in as many data points as particles in the sys- 
tem. In Fig. [6] this data is averaged in bins of size two 
particle diameters. The fluctuations within each bin are 
shown by the vertical error bars. Both curves (for yi and 
y e ) can be relatively well approximated by parabolas: 

Ay( yi ) = a' + b'yi + cyl Ay(y e ) = a + by e + cy\ (1) 

It is obvious that both curves cannot obey exactly the 
parabolic behavior as yi and y e are related by y e (yi) — 
yi + Ayfjji). However, in the cases presented in this sec- 
tion, obtained by slow deposition, the value of Ay is rel- 
atively small compared to yi so that y e (yi) is very close 
to a straight line, leading only to a very small horizon- 
tal shift. This behavior is typical for intermediate Bond 
numbers whereas in the limiting cases no noticeable de- 
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FIG. 6: (Color online) Size of avalanches depending on ver- 
tical position for Bo g — 10 3 . Here the size is measured by 
Ay, the average total downwards motion of a particle after 
deposition (initial position minus final position). On the hor- 
izontal axis the initial position yi (blue squares) and final po- 
sition y e (red circles) are plotted. This leads to two slightly 
shifted curves as y e < y%. The solid lines represent parabolic 
fits Ay( yi ) = -1.8 - 0.31?/; + 0.00035?/? (black, full line) and 
Ay(y e ) = -10.2 - 0.31y e + 0.00038^ (violet, dashed). Addi- 
tionally shown is a fit by Ay(y e ) = —ay e (l — y e /H) predicted 
by the considerations in sec. IVTT1 leading to a « 0.39 green, 
dashed-dotted). 



pendence of Ay on the vertical position could be found. 
For Bo g = a small constant value, below the particle 
diameter (around 1.5 particle radii) is observed. In the 
case Bo g — )• oo no collapse happens, i.e. all Ay = 0. 
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FIG. 7: (Color online) Collapse of the size of the avalanches 
for two different system sizes can be obtained scaling both 
axes by the system height (here: Bo g = 10 3 ). Under the 
assumption of a parabolic profile this scaling leads to a 1/L 
dependence of c (pre- factor of quadratic term in eq. [1]). 

The parabolic behavior can be reproduced also for 
other system heights. In Fig. [71 two different system sizes 
again for Bo g = 10 3 are shown collapsed by scaling both 



axes by the system height H. From this scaling one can 
deduce the system size dependence of the pre-factor of 
the quadratic term in eq. (pQ). The scaling becomes: 

Ay(y, H)=H- f(y/H) cx H ■ (y/H) 2 cx l/H (2) 

when assuming that Ay cx y 2 (parabolic behavior, see 
eq. [Tj) . This l/H dependence could be verified by fitting 
the curves in Fig. [7J Note that the parabolic shape was 
also found when varying the friction coefficient \i and the 
rolling friction coefficient \i r . 
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FIG. 8: (Color online) The histogram of the size of avalanches 
| Ay | for Bo g = 10 3 follows basically a Gaussian. Deviation 
from this behavior can almost fully be suppressed when re- 
moving the bottom and top part of the system. 

Whereas the average of the avalanche size Ay as func- 
tion of the vertical position shows a parabolic profile of 
reasonable quality, there are of course large fluctuations 
around this value. In Fig. [8] we show the distribution of 
the avalanche sizes (here \Ay\) for the entire system, i.e. 
independent on the vertical position. When removing the 
upper and lower part of the system to decrease bound- 
ary effects we obtain a Gaussian distribution, i.e. we get 
an estimate of a typical avalanche size. This typical size 
decreases with increasing Bond number, and in the limit 
of Bo g — )• cx), where no avalanches occur, it vanishes. In 
the limit of Bo g = (no cohesion) the behavior is differ- 
ent, an exponential decay is obtained (Fig. EJ). Here the 
boundaries have no effect, i.e. we get the same behavior 
when removing the upper and lower part of the system as 
done previously. For this Bond number typically the sur- 
face of the structure during deposition grows relatively 
flat, so that large \Ay\ are unlikely as expressed by the 
exponential decay. Due to the monodispersity this sur- 
face is locally almost a flat crystalline surface with heaps 
which consists of a few particles only, in most cases one 
particle. When a particle is deposited on a one particle 
heap it rolls off to rest eventually as a "crystalline" neigh- 
bor aside the particle, resulting in \Ay\ between one and 
two. This leads to the very small range of \Ay\ with con- 
stant probability in Fig. [9] Taking a slightly polydisperse 
system this region would disappear. 
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FIG. 9: (Color online) The histogram of avalanches \Ay\ for 
Bo g — basically shows an exponential decay. Deviation 
from this behavior can be found for \Ay\ between 1 and 2 
particle radii where the probability is about constant. This 
effect cannot be suppressed when removing the bottom and 
top part of the system (as, e.g. for Bo g = 10 3 ). 



given height y^ e (L x width of the two dimensional system 
in units of particle radii): 

N dJ (y i , e )=L x dy'p dJ (y') (3) 
Jo 

The final position y e of particles can be related to the 
position yi of deposition by the avalanche profile Ay: 

ViiVe) =Ve- Ay(ye), or y e {y%) = Vi + ^y(yi) (4) 

In this notation Ay is negative as the motion of the par- 
ticles is downwards (due to gravity). Therefore, yi is 
larger than or equal to y e . As particles are never de- 
stroyed the number of particles deposited up to a given 
height, Nd(yi), will stay the same, but shifted to a lower 
height, Nf(y e ), where yi and y e are related by eq. (|4]). 
Together with eq. (|3]) this leads to: 

N f (y e )/L x = N d (y t (y e ))/L x = f ^ dy' p d (y') (5) 

Jo 

=G(yi(y e )) 



In this section we studied the collapse of the struc- 
tures occurring in small avalanches we analyzed statis- 
tically. We suggest to characterize these avalanches by 
their "size" , showing a typical dependence on vertical po- 
sition, a parabolic shape for the specific systems investi- 
gated in this section. In the following section we will use 
this characteristic behavior to be able to derive the final 
density profile from the "deposition density" . In section 
PVTl the same concept will be shown to be applicable also 
for other protocols of generating loose structures. 



V. THEORETICAL ANALYSIS OF THE 
AVALANCHES 

In the previous sections we mentioned that the dy- 
namics leading to a final configuration is determined by 
small avalanches occurring during the deposition process. 
All these compaction events contained in the function 
Ay(y e ), which is given by the difference between the ini- 
tial position yi and final position y e . Note that Ay can be 
plotted (e.g. fig. [6]) as function of the final position y e or 
alternatively as function of the position of deposition y^. 
The aim of this section is to relate the final density profile 
to the dynamic process of deposition and collapse by us- 
ing Ay(y e ), showing how the avalanches produce the final 
density pf(y e ) from the deposition density p d {yi). The 
deposition density is defined by the number of particles 
deposited within a volume. As the structure collapses 
between the depositions the deposition density is not in- 
dependent on the collapsing, and it is possible that at (al- 
most) the same position several particles are deposited. 
Thus, locally within a fixed volume even more particles 
could be deposited than typical for a dense packing. 

We first calculate the number N d j of particles up to a 



This relates Nf to the deposition density whereas eq. (j3j) 
relates Nf to the final density. The function G here is 
formally introduced for the integral as abbreviation, by 
derivation of G the density is retrieved. The final density 
can be obtained by derivation of Nf / L x using eq. ([5]) : 



p (y )= - N f( ye " j d C(y (y )) dG ^ dyi 
Pf Ve dy e L x dy e V% Ve dyi dy e 



= Pd(yi(ye))^- 
dy e 

= Pd{vi{ve)) fi- 



dAy(y e ) 
dy e 



The deposition density p d (yi(y e )) in principle can be ex- 
pressed directly by y e introducing p f d (y e )- As usually the 
functional behavior of both functions is not known, but 
only values for specific yi and y e , the transformation can 
be done for each point by simply using eq. ((4]), i.e. re- 
placing each yi by y e = yi + Ay{yi). Summarizing, to 
calculate the final density profile one needs to know the 
deposition density p d and the avalanche profile Ay. Note 
that the avalanche profile dependence on both y e and yi 
is needed, which can be calculated from each other for 
some cases as shown later. For experimental situations 
these quantities are not known. However, the relation 
between pf and p d (eq. [6]) can be used to calculate the 
deposition density from the final density in the slow de- 
position limit, when assuming a parabolic profile as found 
in the simulations before. 

In figure [10] we use eq. (j6j) to calculate the final den- 
sity from the deposition density by using the parabolic 
fit for Ay (fig. [6]). In practice first the deposition den- 
sity curve is shifted on the horizontal axis by y e = 
yi — {a ,J rb , yi J rC , yj) 1 then multiplying the deposition den- 
sity with the right hand side of eq. (|6)), 1 — (6 + 2q/ e ), i.e. 
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FIG. 10: (Color online) Using the parabolic approximation 
(fig. E]) for the average avalanches the final density (here vol- 
ume fraction vf) can be calculated from the deposition den- 
sity (here shown: volume fraction Vd for Bo g = 10 3 , inset). 
There are strong fluctuations in the deposition density which 
are induced by the irregularity of the avalanches. To obtain 
a smooth curve we use a fit function (here: power law fit, 
with exponent 0.15) to calculate the final density (black line) 
matching relatively well with the measured density profile for 
sufficiently large y (except close to the bottom) 



using the derivative of Ay(y e ) which is a linear function. 
If the deposition density would be constant this would 
lead to a linear profile for the final density. However, the 
deposition density is not constant, explaining the non- 
linear behavior for the final density. Additionally the 
deposition density shows strong fluctuations, but by as- 
suming the avalanches to follow the averaged parabolic 
behavior the corresponding fluctuations in the avalanche 
profile are not included. The calculated curve matches 
relatively well the profile measured in the simulations for 
sufficiently large values of the vertical position. Close to 
the bottom, however, the calculated curve deviates from 
the measured one. In this region the deposition density is 
very small, i.e. almost the one of pure ballistic deposition. 
This can be understood as the system needs to gain a suf- 
ficient amount of weight for the collapse to start (cf. also 
sec. IVip . This should correspond to a higher initial slope 
of Ay(y e ) which is not reflected in the parabolic approxi- 
mation (eq.HJ. In this region a higher order terms would 
be necessary for reproducing also the system bottom. 

The same analysis has been done also, e.g. for Bo g = 
10 2 as shown in Figs. [11] and [T2l In this case the de- 
position density shows somewhat lower fluctuations as 
for Bo g = 10 3 (Fig. fTUj) . To quantify this we esti- 
mated the fluctuations of the deposition density at ver- 
tical position y = 200 for both cases. For Bo g = 10 2 
we obtained around 15% whereas we estimated around 
20% for Bo g = 10 3 . For the case of Bo g — )• oo there is 
no avalanching at all (cf. sec. IIV|h and trivially the fi- 
nal density equals the deposition density. This is very 
similar for very large Bo g , but as some avalanches occur 
there are some relatively small fluctuations in the aver- 



FIG. 11: (Color online) Size of avalanches depending on ver- 
tical position for Bo g — 10 2 . Here the size is measured by 
Ay, the total downwards motion of the particle after deposi- 
tion (initial position minus final position). On the horizontal 
axis the initial position yi (blue squares) and final position 
y e (red circles) are plotted. The lines represent the parabolic 
fits Ay( yi ) = -1.7 - 0.24^ + 0.00056?/ 2 (dashed, violet) and 
Ay(y e ) = -4.7 - 0.23y e + 0.0005% 2 (full line, black). 
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FIG. 12: (Color online) Using the parabolic approximation 
(fig. Hip for the average avalanches the final density can be 
calculated from the deposition density (here shown: volume 
fraction for Bo g = 10 2 ). There are strong fluctuations in the 
deposition density which are induced by the irregularity of the 
avalanches. As before (for Bo g — 10 3 ), to obtain a smooth 
curve we use a fit function (here: power law fit, with expo- 
nent 0.13) to calculate the final density (black line) matching 
relatively well for sufficiently large y (except close to the bot- 
tom). 



age profile. Away from this limit, but still close enough 
that the density profile is very similar to the Bo g — >> oo 
case, as e.g. for Bo g = 10 , very large fluctuations in the 
avalanche profile are observed. Thus, the parabolic pro- 
file cannot easily be identified. Still the theory works well 
as the final density is very close to the deposition den- 
sity so that even a very inaccurate fit for the avalanche 
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profile does not affect the calculated density profile too 
much. In the limit of Bo g = the avalanche profile is 
a constant (cf. sec. HV]h i.e. all grains are slightly shifted 
downwards by the same amount (except boundary effect 
at the bottom). As then the derivative vanishes the final 
density equals the deposition density. 

Here, we showed how the parabolic avalanche profile 
can be used to calculate the final density profiles from the 
deposition density in the case where gravity acts during 
deposition. In the next section the same concept will be 
used to the simpler case of collapse of the system after 
deposition is complete. These two cases could be then 
related in section IVTTI 



VI. COLLAPSE AFTER DEPOSITION 
COMPLETE 

In the previous sections we investigated the case where 
gravity acts during deposition, leading to relatively com- 
plex shape of the density profiles and a parabolic char- 
acteristics of the avalanche size. For this case we showed 
that these avalanche profiles can be used to relate the fi- 
nal density profile to the deposition density. In this sec- 
tion we will analyze the case when the particles are first 
deposited, then gravity is switched on and the structures 
collapse. This case is even simpler and can later be used 
to understand the more complex system studied before. 
In this case the initial density pi characterizes the sys- 
tem (instead of the deposition density as in the previously 
discussed situation) . Using the off lattice version of bal- 
listic deposition as presented in Refs. [4l|, |42| with stick- 
ing probability one, vertically falling particles stick when 
they touch an already deposited particle. This leads to 
a fixed initial density. Lower densities can be obtained 
by using a capture radius r capt i.e. particles stick to each 
other when they are within a certain distance during the 
falling of the depositing particle. More precisely: When 
the distance between the center of masses of two particles 
is below 2 • r capt the particles stick and the falling particle 
is pulled along the connecting line towards the already 
deposited particle. This capture radius is a measure for 
the distance between the branches of the deposit and the 
resulting density is inversely proportional to r cap t [HI, 
^capt = 1 gives the original method. The resulting ini- 
tial structures are shown in Fig. [13] These structures 
obtained with different capture radius will be later used 
to study the influence of the initial density. 

First we will investigate the behavior using r cap t = 1. 
Figure [14] shows the density profile before the collapse 
which is the same that we got in the limit of Bo g — >> oo 
in sec. IIII1 also independent on vertical position. After 
this deposition is complete gravity is "switched on" and 
the structure abruptly collapses. Here we choose a Bond 
number of Bo g — 10 3 . This leads to a final structure 
with higher density, in this case also independent on the 
vertical position (Fig.fT4|). As no particles are added after 
the initial deposition the final system height is lower. 



Similarly as we did before we analyze the size of the 
avalanches Ay as defined in sec. HVi Figure [15] shows 
a linear dependence of Ay on either y e and yi. The fit 
parameters of the two lines can be related to each other 
by the relation between yi and y e (eq. [4]). Assuming 
Ay(Ve) = a — by e and Ayfjji) = a' — b'yi the values a' 
and b' can be calculated from a and b (see app. [AJ as: 



The vertical dependence of Ay can be used similarly as 
before to calculate the final density from the initial den- 
sity by using eq. (|6]). The calculated density profile us- 
ing this linear dependence reproduces the obtained final 
density profile very well as shown in Fig. [TU In this case 
the agreement is better as now the initial density is not 
fluctuating very much in contrast to the cases discussed 
in sec. [V] The density increase Ap (or volume fraction 
increase Av) can be directly calculated by the constant 
slope of Ay(y e ): 

Ap Av dAy(y e ) 

— — (8 J 

Pi Vi dy e 

For the same parameters (Bo g — 10 3 ) we studied the 
effect of the system height H on the density increase 
while still keeping the initial density fixed (fig. [T6]) . A 
logarithmic fit matches the data best. This fit certainly 
cannot continue to infinity as there is a limit for the den- 
sity p max given by the random close packing (see also fig. 
E]), leading to a (Az//V;)max of 1.19(~ /WxMni - 1). 

Using initial capture radii as described above we study 
the influence of the initial density on the relative density 
increase Avjvi (fig. |TTJ) . We could obtain the best fit 
when using a power law with exponent of about 1.64. 

We showed in this section that the linear avalanche 
profile is a characteristic feature of compacting from a 
depth independent to a depth independent structure, ob- 
tained here for systems generated by ballistic deposition 
collapsing due to gravity. More complex avalanche pro- 
files with non-constant derivative will transform homoge- 
neous structures into inhomogeneous structures. Thus, 
we expect the linear profile to be obtained in all cases 
where a homogeneous initial system compacts to a ho- 
mogeneous final system. These homogeneous compaction 
processes are investigated in different research areas as 
e.g. discussed in Refs. (53l-[57j. In addition in the next 
section we will show that also for the more complex pro- 
cess when gravity acts during deposition (sec. IIII|) this 
linear profile can be used to derive the parabolic profile 
of the avalanches. 



VII. RELATION BETWEEN DEPOSITION 
UNDER GRAVITY AND SWITCHING ON 
GRAVITY AFTER DEPOSITION 

For the very fast process a linear profile for Ay de- 
pending on vertical position has been found (cf. fig. [T5j) 




FIG. 13: Initial structures generated by ballistic deposition with increasing capture radius r ca pt. 




FIG. 14: (Color online) The initial and final density are about 
constant when depositing first and then collapsing the system 
(Bo g = 10 3 ). Using the linear dependence of the avalanches 
Ay on the vertical position (see Fig. I15[) the final density can 
be calculated from the initial density using eqs. (4|) and (|6]). 
The results support the analytical considerations. 



FIG. 15: (Color online) The linear dependence of the 
avalanche sizes Ay(yi, e ) explains the homogeneous density 
increase as seen in fig. 1141 The linear fits are Ay(y e ) = 
2.2 - 0.62y e and Ay( yi ) = 1.9 - 0.3%. 



whereas the slow deposition limit shows a parabolic 



profile for Ay depending on vertical position (cf. figs. 
I6|lip . In this section we will discuss how a relation be- 
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FIG. 16: Dependence of volume fraction increase Avjvi on 
system height H when first deposited and then collapsed. A 
logarithmic fit matches the data best (here y = —0.94 + 
0.26 ln(x)). The limit of random close packing defines the 
largest possible value for Av/vi of 1.19, which will be ap- 
proached for infinite system heights. 
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FIG. 18: Sketch illustrating the procedure of depositing the 
grains slice by slice. The first slice deposited is compacted 
by internal collapse. The same is true for each "freshly" de- 
posited slice. The slices below are compacted by the added 
weight of the slices above. Periodic boundary conditions in 
horizontal direction are imposed (illustrated by dashed lines). 
The figure also illustrates the definition of the symbols used 
here. The n slices are numbered from 1 to n. A slice i col- 
lapses from pi-i to pi while its height decreases from hi-i to 
hi, where hi = pi-i/pihi-i. 



FIG. 17: Dependence of volume fraction increase Avjvi on 
volume fraction Vi of the initial system (system first deposited 
and then collapsed). Different densities could be reached by 
increasing the capture radius for ballistic deposition (fig. I13p . 
A power law fit with exponent 1.64 fits relatively well (power 
law fit results in y = 0.158x" 164 ). 



tween both can be established. By this relation also the 
parabolic profile is put onto a more fundamental basis 
like the linear profile for the homogeneous collapse. 

Let us imagine depositing particles slice by slice as 
sketched in Fig. [181 The slices are thin parts of the sys- 
tem in vertical direction spanning the full system width in 
horizontal direction. They can be considered as systems 
with very small initial height ho. In each slice the depo- 
sition will be immediately followed by the collapse. How- 
ever, there will be not only an "internal collapse" within 
the "freshly" deposited slice, but also a compaction of the 
slices below due to the additional weight of the "freshly" 
deposited slice. 



Let us first consider systems composed of a small num- 
ber n of slices. The case n = 1 (one slice) is the same as 
discussed in the previous section: The system collapses 
"internally" leading to an increase of the density from po 
to pi while the height reduces from ho to h\. Here we de- 
note the slice number as 1 (cf. Fig. [15]). As shown in the 
previous section the avalanche sizes have a linear profile 
Ay(yi 1>} ) — Siyi r \ Si is the slope in slice 1, and is the 
same for all freshly deposited slices when the height ho 
is kept constant. The vertical position within slice 
1 is measured from its bottom (y^ = 0.../ii). This 
notation will be used in the following for each slice i: 
yf* = . . . hi. The case n = 2 (two slices) means adding 
an additional slice to the case n = 1. Then the lower 
slice (slice 2) experiences an additional compaction by the 

added weight expressed by the corresponding avalanche 

(2) 

size C2ye assuming a linear behavior for this relatively 
fast process similar as for the internal collapse. This is 
justified at least for the limit of small slices later consid- 
ered in this section. The upper slice (slice 1) will be com- 
pacted internally and additionally will move downwards 
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by C 2 h 2 (= hi — h 2 ) as the slice below is compacted. 
Summarizing for the two slices we get: 



AyW(yW) = S lV ^ + (9) 

hi —hi 

Ay^(yW) = C 2 y^ + S 1 (l + C 2 )y^ 



For slice 2 the internal compaction from the first step 
S\y^ has been transformed by using that p 2 = (1 + 
C 2 )pi (cf. eq. [8]), leading to h 2 = pi/ ' p 2 h\ = 1/(1 + 
C 2 )hi. Adding a further slice leads to the case n = 3, 
where the two slices are compacted due to the additional 
weight. Each of these compactions is accompanied by a 
downwards shift of the slices above. This leads to: 



Ay^\y^) 
Ay^(y^) 



Siy™ + C 2 h 2 + C 3 h 3 (10) 

S v ' 

(hi—h 2 ) + (h 2 —h3)=hi—h3 

C 2 yi 2) + 3,(1 + C 2 )yP + C 2 h 2 +C 3 h 3 
from step 2 

Csy^ 

+ C 2 (l + C 3 )y^ + + C 2 )(l + C s )yi 3) 

V v ' 

from step 2 



Imagining continuing this iterative procedure, one ob- 
tains the case of n slices. For the top slice this results 



A^ (1) (^ 1} ) = S iy ^ + C 2 h 2 + C 3 h 3 + • • • + CM U) 

hi—hn 



The first term is the internal collapse where the other 
terms are the shift due to the compaction of all slices 
below (2 to n) in this last step. For the bottom slice we 
get: 



AyW(yM) = 5i(l + C 2 )(l + C 3 )-...-(l + C n )j/W 

C n yW+C n - 1 (l + C n )yW + ... (12) 
+C 2 (1 + C n )(l + C n _i) ....-(! + C 3 )yi n) 



Here all terms represent a collapse in the slice either inter- 
nally by its own weight when deposited in the first step, 
or when collapsing due to added weight in the following 
steps. In addition these collapses have to be transformed 
to a yi n ^ dependence (see above). For an arbitrary slice 
i somewhere in the system we get both types of terms as 



in eqs. §TQ) and (fT2j) : 

Ay«(y«) = S 1 {l + C 2 ).....(l + C i )y® 

CiyW+Ci-iCl + COi/W (13) 
+ --- + C 2 (l + C 3 )-...-(l + C,)^ ) 

-\-C n h n + C n -\h n -i + • • • + Cn-i+ihn-i+i 

V v ' 

h n -i-h n 

-\-C n —\h n —i + ■ ■ • + C n -ih n -i 



-i-h n 



■C 2 fc 



2^2 



hi—hi+i 



The part of the expression independent on yi^ represents 
the shift due to compaction by the weight of the above 
added slices in n — i steps considering all slices below. It 
consists of n — i times i — 1 terms and can be written 
shortly as: 



'shift 



(14) 



The limit of large n while keeping the total system 
height constant gives very small slices where the part 

(i) 

A?/g hift dominates as for very small systems the internal 
collapse almost vanishes (cf. fig. [T6|) . Therefore, in the 
following we will only consider this term to show that 
we approximately obtain a parabolic behavior. Let us 
assume that the hj are linear in i: 



a- I h , 

n 



a < 1 



(15) 



This means that deeper in the system (larger i) the width 
of the slice is smaller. Note that for the case a«l this 
can be understood as a linearization. This case means 
that the overall compaction is not large as it is the case 
for intermediate Bond numbers. From eqs. ([T5]) and (fT4|) 
we obtain: 



3=1 



J , 

1 — a ) ho 



a— ho = (n — i)a — ho 
n n 



(16) 



This is a quadratic dependence on the slice number i. To 
compare to our results we have to transform i to verti- 
cal position y e , which is obtained when summing up the 
height hi of all slices: 



i+l 



'3+i 



(17) 
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Using the approximation (fT5]l we obtain: 



5 (i) ~ Wl- fl 



(18) 



= fto [n- ^(n 2 -n) -i/i [l-a(n-l/2)](19) 

The detailed derivation is given in appendix [B] From 
this equation we can obtain i(y e ): 



n — a/2(n 2 - n) 
h [l-a(ra-l/2)] + 1 - a{n - 1/2) 



(20) 



We assume that we are in the limit of relatively small a. 
Neglecting all terms in a in eq. ([20]) corresponds to ne- 
glecting terms in a 2 in eq. JT5J). With this simplification 
additionally using ho = H/n we obtain i = n — y e /ho = 
n(l — y e /H), leading to: 



A2/ S hift (2/e(*)) 



[™ - i(ye)]i(y e )^ 



ay e i 



2/e\ 



(21) 
(22) 



This behavior is plotted in fig. [6] (green curve). Note that 
in this figure the Ay is negatively defined as opposed to 
the definition used in this section. This curve fits rela- 
tively well the measured curves except coming close to 
the top. This can be explained by the existence of a 
small "crust", i.e. an accumulation of particles at the 
top of the system in the simulations which is not con- 
sidered in the analysis in this section. Probably this is 
also the reason for the slightly different pre-factors of the 
parabola: From eq. ([22]) we obtain a = 0.39 leading to a 
pre- factor of the quadratic term of a/H = 0.00049 which 
is somewhat larger than the value obtained previously of 
0.00038. The value of a = 0.39 is at least reasonably 
small to ensure that the considerations of this section 
agree roughly with the simulation results. Previously by 
scaling Ay for different system sizes we obtained that 
the pre- factor of the parabola scales as 1/H (cf. eq. [2]). 
This implies that a is independent on system size H, for 
each specific Bond number, additionally indicating by its 
value how good the approximations of this section are. 

Thus, in the limit of small a we could show that the 
linear behavior of Ay when collapsing after deposition 
complete leads to a parabolic behavior when collapsing 
during deposition. Note that this a represents the differ- 
ences in heights of the top and the bottom slice, i.e. the 
assumption of small a is true when the density difference 
between the density close to the bottom and at the top 
is small which is the case in all our cases studied here (cf. 
Fig. [2]). In the structures studied within our model in the 
previous sections (see e.g. sec. IIII|) a small "crust" (par- 
ticle accumulation) at the top leads to a density increase 
again. This will lead to a shift in the parabolic profile 
to the right (to the top). As we discussed previously the 
deposition/collapse process is not continuous, so that the 
parabola is only an average of a very noisy distribution 



of Ay. Additionally, the deposition density is not con- 
stant, but slightly increasing (cf. figs. [TUl and [T2|) accom- 
panied by relatively large fluctuations. For these reasons 
we can only expect a rough matching of our theory with 
the simulations. Nevertheless the parabolic behavior has 
been observed relatively clearly. 



VIII. CONCLUSION/OUTLOOK 

We studied the generation of fragile granular struc- 
tures by a deposition/collapse process. In one extreme 
case where the deposition is sufficiently slow to allow the 
system to collapse and relax due to gravity after the depo- 
sition of each single grain we studied the influence of the 
granular Bond number on the density profile. For inter- 
mediate Bond numbers the density decreases with height 
due to the compaction of the powder's own weight. We 
studied the generation process dynamics which is discon- 
tinuous in small avalanches. These avalanches showed a 
parabolic behavior and can be used to calculate the fi- 
nal density profile from the deposition density. In the 
other extreme case of collapse after deposition complete 
we found that the density is constant with vertical po- 
sition, and that the avalanche size depends linearly on 
vertical position. We could relate the parabolic behav- 
ior to the linear one by imagining a slice by slice de- 
position/collapse process. Note that the linear behavior 
investigated here for the case of ballistic deposition fol- 
lowed by a gravitational collapse will be found for all col- 
lapse/compaction processes of homogeneous initial struc- 
tures to homogeneous final structures. Therefore the 
concept of avalanches introduced in this paper is of gen- 
eral applicability to granular structures collapsing due to 
gravity or similar forces. 

Our results maybe directly verified experimentally, as 
already mentioned in the introduction, e.g. by using a 
Hele Shaw cell [43-45] which can be tilted to effectively 
change gravity. To apply the model presented here more 
specifically, e.g. for snow compaction, more realistic mi- 
croscopic properties including aging processes would have 
to be used. For cake formation processes, instead of 
gravity a porosity dependent drag force could be ap- 
plied. In this context an explicit consideration of the 
pore fluid/gas could be needed. The influence of the pore 
fluid/gas should be in particular studied for the fast com- 
paction process presented in this paper. 
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Appendix A: Relation between slopes 

The linear dependence of avalanches is found as well 
in y e as in yi (see fig. [T5|) . In this section the relation 
between the two lines is derived in detail. Assuming 

Ay(y e ) = a-by e and Ayfa) = a' - b'yi (Al) 

the values a' and b r can be calculated from a and b as 
shown in the following. The relation between yi and y e 
can be written as: 



yi = y e - Ay(y e ) = 



y e -(a- by e ) 
y e (l + b) - a 

1+6-6 

a 1 



1 + 6 



1 + 6 
b 



Vi 



1 



1 



(A2) 
(A3) 

(A4) 
(A5) 



According to © and (|A1J| y e can be written as: 

Ve = Vi + (a' ~ b'yi) (A6) 
Comparing (|A5[) and ()A6[1 results in: 
6 , a 



b' 



1 + 6' 



1 



(A7) 



From this or by a similar derivation the inverse relations 
can also be obtained: 

(A8) 



l-V 



l-V 



Appendix B: Derivation of y e (i) in linear 
approximation for hi 

Here we show the details of the derivation to obtain 
eq. (jUJ) from eq. (fT8|) : 



Ve(i) - Yl ( 1 ~ a 

3 = 1 



j + i 



h 



( 



ha 

= h {n - i) 

n 



\ 



n—i 



V(n-i)t (n-t+l)(n-t)/2/ 

= fto(^ — — hoa(n — i)i/n — hoa(n — i + l)(n — i)/2 



h (n - i) - h a -(n + 1) - i ((n + l)/2 + n/2) 



n — t^ 77,2 — n ) — [1 — a ( n — 1/2)] 
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